function F = solve_eqm_cf2_outer(param, const, eqm)
%% Counterfactual: Export Promotion w/ Balanced Trade  

    % Initial quality
    q_indicator = eqm.q_indicator;
    
    % Initial integrals  
    integ = eqm.integ;
    
    
%% Solve the open economy model with quality choice

    % Setting
    tol = 1e-12;
    maxIter = 100;
    maxHistory = 5;
    maxRepeat = maxHistory+1;
    
    % Placeholder
    q_opt_index_history = zeros(param.Nf,maxHistory); 
    switch_id_history = zeros(param.Nf,maxHistory+1);
    switch_qindex_history = zeros(param.Nf,maxHistory+1);  
    
    % Initialize
    eqm_cf.Pi0 = ones(1,param.Q);
    eqm_cf.Pi1 = ones(1,param.Q);  
    nq_new = cal_nq(const, q_indicator);
    fix_flag = 0;
    repeat = zeros(1,maxHistory);
    
    % Figure
    figure; 
            
    % Iteration
    err = 1;
    iter = 0;
    while (err >= tol && iter < maxIter)
        
        % save previous H quality choice
        if iter > 1
            q_opt_index_history = [q_opt_index, q_opt_index_history(:,[1:end-1])];
        end   
        
        % save previous H switch id/qindex
        if max(repeat(2:end)) > 0           
            temp_switch_id = zeros(param.Nf,1);
            temp_switch_id(switch_id) = switch_id;
            switch_id_history = [temp_switch_id, switch_id_history(:,[1:end-1])];            
            temp_switch_qindex = zeros(param.Nf,1);
            temp_switch_qindex(switch_id) = switch_qindex;
            switch_qindex_history = [temp_switch_qindex, switch_qindex_history(:,[1:end-1])];
        end
        
        % new guess
        Pi0 = eqm_cf.Pi0;
        Pi1 = eqm_cf.Pi1;          
        nq = nq_new;
        iter = iter + 1;
                
        % solve equilibrium
        eqm_cf = solve_eqm_cf2_inner(param, const, integ, eqm); 
                          
        % choose new quality 
        [q_opt_index, ~] = grid_search(param, const, eqm_cf.Pi0, eqm_cf.Pi1, eqm_cf.ExportCutoffQ, eqm_cf.ExportProbQ);        

        % fix quality after maxRepeat
        if max(repeat(2:end)) >= maxRepeat 
            if fix_flag == 0
                % find id and qindex
                fix_id = nonzeros(max(switch_id_history,[],2));
                fix_qindex_history = switch_qindex_history(fix_id,:);
                fix_qindex = max(fix_qindex_history,[],2); 
                % set flag
                fix_flag = 1;
                repeat = repeat*0;
            elseif fix_flag > 0
                % find id and qindex
                fix_id2 = nonzeros(max(switch_id_history,[],2));
                fix_qindex_history2 = switch_qindex_history(fix_id2,:);
                fix_qindex2 = max(fix_qindex_history2,[],2);  
                % append
                fix_id = [fix_id; fix_id2];
                fix_qindex = [fix_qindex; fix_qindex2];
                fix_qindex_history = [fix_qindex_history; fix_qindex_history2];
                % set flag
                fix_flag = fix_flag + 1;
                repeat = repeat*0;
            end
        end
        if fix_flag > 0
            q_opt_index(fix_id) = fix_qindex;                
        end           
        
        % update qindicator
        q_indicator = cal_indicator(param, q_opt_index); 

        % update integrals
        integ = cal_integral(const, q_indicator, eqm_cf.ExportProbQ); 
        
        % update n(q)
        nq_new = cal_nq(const, q_indicator);                     
        
        % check convergence 
        err1 = max(abs(nq_new-nq));
        err2 = max(abs(eqm_cf.Pi0-Pi0));
        err3 = max(abs(eqm_cf.Pi1-Pi1));        
        err = max([err1, err2, err3]);
        
            % visualize progress
            scatter(const.Qgrid, nq_new, 'Marker', '.')
            xlabel('q')  
            title(strcat('iter=',num2str(iter)))
            drawnow limitrate
        
        % firm types changing quality
        switch_id = find(q_opt_index~=q_opt_index_history(:,1));
        switch_qindex = q_opt_index(switch_id);
        switch_num = size(switch_id,1);
        
        % repeat history
        for i = 1:maxHistory
            repeat(i) = repeat(i)*min(q_opt_index == q_opt_index_history(:,i)) + min(q_opt_index == q_opt_index_history(:,i));
        end                                  
        
            % print
            fprintf('iter%.0f: %.0f/%.0f types change quality', iter, switch_num, param.Nf); 
            if switch_num > 0
                if max(repeat(2:end)) > 0
                    fprintf(', repeat ='); fprintf(' %.0f', repeat); 
                end
                if switch_num <= 5
                    fprintf(', id ='); fprintf(' %.0f', switch_id);
                    fprintf(', qindex ='); fprintf(' %.0f', switch_qindex);
                end 
            elseif switch_num == 0
                fprintf(', err = %e', err);
            end
            fprintf('\n')            
            
    end
    close
    
    
    
%% Return

    eqm_cf.iter_outer   = iter;
    eqm_cf.err_outer    = err;
    if fix_flag > 0    
        eqm_cf.fix_id               = fix_id;
        eqm_cf.fix_qindex           = fix_qindex; 
        eqm_cf.fix_qindex_history   = fix_qindex_history;
    else
        eqm_cf.fix_id               = nan;
        eqm_cf.fix_qindex           = nan;         
        eqm_cf.fix_qindex_history   = nan;
    end 
    
    eqm_cf.q_index     = q_opt_index;
    eqm_cf.q_indicator = q_indicator;
    eqm_cf.integ       = integ;
    
    eqm_cf.nq               = nq_new;
    eqm_cf.first_nonzero    = find(eqm_cf.nq, 1, 'first');
    eqm_cf.last_nonzero     = find(eqm_cf.nq, 1, 'last');  
    eqm_cf.first_nonzero_q  = const.Qgrid(eqm_cf.first_nonzero);
    eqm_cf.last_nonzero_q   = const.Qgrid(eqm_cf.last_nonzero); 
    
    eqm_cf.density_all      = const.density_omega;  
    eqm_cf.density_exp      = const.density_omega.*full(sum(eqm_cf.ExportProbQ.*eqm_cf.q_indicator,2));
    eqm_cf.density_nonexp   = const.density_omega - eqm_cf.density_exp;   
    
    eqm_cf.density_exp_ctn = min(eqm.density_exp, eqm_cf.density_exp);
    eqm_cf.density_nonexp_ctn = min(eqm.density_nonexp, eqm_cf.density_nonexp);
    eqm_cf.density_switcher = eqm_cf.density_all - eqm_cf.density_exp_ctn - eqm_cf.density_nonexp_ctn;    
    
    F = eqm_cf;

    
end